Dynamical behavior of a lattice glass model on a random graph: comparison with 

Mode Coupling Theory 
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We study the dynamical behavior of a lattice model of glass former on a random graph, where 
no corrections to the mean field description are expected. We find that the behavior of dynamical 
correlation functions and dynamical susceptibility are consistent with the quantitative predictions 
of the Mode Coupling Theory of the glass transition. 
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The glass transition is characterized by a dramatic rise 
of the relaxation times of the system, over a small interval 
of temperatures or densities. It has long been supposed 
that this phenomenon is connected to the growing of a co- 
operativity length in the dynamics [1, 2]. This idea sug- 
gests the existence of dynamical heterogeneities in super- 
cooled liquids, that is groups of particles whose motion 
is spatially correlated, on a scale that becomes large as 
the glass transition is approached. Such heterogeneities 
were identified in experiments [3] and simulations [4]. A 
way to quantify this dynamical heterogeneities is in terms 
of the four point (dynamical) susceptibility X4(0 [5, 6], 
defined as Xi{t) = V [{(t){tf) - (0(0)^], where 0(0 = 
p(r, 0)p(r, t) — (p)^, and p{r, t) is the density at position 
r and time t. The susceptibility X4(0 ^he volume in- 
tegral of a four point correlator, Xi{^) = j G4,{r,t)(fir, 
where 



G4(r,t)=(p(0,0)p(0,t)p(r,0)p(r,t)} 

-(p(0,0)p(0,t))(p(r,0)p(r,t)}. 



(1) 



measures the correlations in space of local time correla- 
tion functions. The function X4(t) has been computed 
in many model glass formers, and indeed it displays a 
maximum at a time t* ~ Tq, where Tq is the relaxation 
time of density fiuctuations, with the maximum X4(^*) 
growing when approaching the transition, signaling the 
growing of the range of G4(r,t*), that is of dynamical 
correlations. 

A mean field level description of the glass transition 
is the Mode Couphng Theory (MCT) [7], which makes 
many predictions on the form of the relaxation near the 
transition. In particular it predicts a power law diver- 
gence, Ta ^ \T — Tc\ \ of the relaxation times. It was 
recently shown [8] that, within the mean field descrip- 
tion of MCT, the time t* where X4(^) has a maximum 
diverges with the same exponent of the structural relax- 
ation time, t* ^ Ta ~ \T — Tc\~'^ ■ Furthermore, the value 
of the maximum diverges as X4(^*) ~ t*^, with x = I/7. 
One also finds that, in the early and late /3 regime, the 
dynamical susceptibility behave as X4(0 ^ ^''j respec- 



tively with /i = a and p — b [9], where a and b are the 
standard MCT exponents that describe the approach and 
departure from the plateau. 

To check these results using molecular dynamics on a 
realistic model of glass forming liquid, like a Lennard- 
Jones mixture, is not an easy task for two reasons. The 
first is that simulations are computationally heavy, and 
it is difficult to reach very long simulation times. The 
second is that the predicted Mode Coupling transition is 
smeared out by corrections to mean field, the so called 
hopping processes, so that near the transition a crossover 
to a different regime is observed. A class of models 
in which the MCT is exact are some spin glass mod- 
els with a discontinous transition, for example the p-spin 
model [10]. However, such models are characterized by 
quenched disorder, that gives rise to strong finite size 
effects, due to sample to sample fluctuations of the crit- 
ical temperature, as pointed out in a recent paper [11]. 
For this reason, in such models the asymptotic critical 
behaviour, that is expected only in a very small region 
around Tc, is not easily observed. 

In this paper, we study the dynamical behavior of a lat- 
tice model of glass forming liquids on the random graph 
(or Bethe lattice). Being defined on lattice, the model is 
much less computationally expensive than more realistic 
models. Furthermore, due to the tree-like nature of the 
graph, we expect the critical behaviour to be described 
correctly by Mode Coupling Theory. Finally, the model 
has no quenched disorder (apart from the randomness 
of the graph), so that it is reasonable to expect not too 
large sample to sample fluctuations. The model is de- 
flned as follows [12]. The space is partitioned in regular 
cells, such that not more than one particle can have its 
center of mass inside the cell. The position of the particle 
inside the cell is discretized, so that it can assume a flnite 
number q of positions. The model is therefore described 
by the following Hamiltonian: 



n 



ninj(f)ij{(Ji,(Tj) 



(2) 
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where = 0, 1 represents whether the i-th cell is oc- 
cupied or not, (Tj = 1 ... (7 represents the position of the 
particle inside the cell, and (j>ij {ai,(7j) the interaction be- 
tween two particles in neighboring cells i and j, with po- 
sitions fTj and aj. 




FIG. 1: a) The model in two dimensions: the space is par- 
titioned in square cells, and each cell can be occupied by at 
most one particle in one of four positions (little circles). A 
particle iu a given position (big shaded circle) forbids the pres- 
ence of another particle in the positions colored in black, b) 
The model on the random graph with 2 = 4. 

A particularly simple realization of the model is shown 
in Fig. la for two dimensions. In this case the space is 
partitioned in square cells, with 9 = 4 positions per cell, 
and the interaction (pij{(Ti, Uj) is zero, if the positions 
and <jj are "compatible" , infinite otherwise. In Fig. lb 
the same model is defined on the random graph, with 
the same value 2; = 4. The graph in this case is extracted 
randomly from the set of all regular graphs where each 
site has exactly z neighbors [13]. This kind of graph is 
locally tree-like, as there are no cycles up to a distance 
scaling as logA^ from a typical site. However, the pres- 
ence of long cycles induces frustration and ensures that 
the system is statistically homogeneous. The mean field 
approach is exact on such graphs, due to the fact that 
the local fields acting on a site are uncorrelated in the 
thermodynamic limit. 

In Refs. [12] the thermodynamics of the model on the 
random graph was studied analytically. Using a 1-step 
replica symmetry breaking (RSB) ansatz, it was found 
[12] that in the case z = 6 the model has a "dynamical 
transition" at a density pd = 0.808, where self-consistent 
1-RSB solutions appear, signaling the appearing of an 
extensive number of metastable states. The transition 
is called dynamical, because it is usually connected to a 
MCT-like dynamical freezing, as it happens in the p-spin 
model [10]. At a higher density pK = 0.822, the free 
energy of the 1-RSB solution equals the free energy of 
the replica symmetric solution, signaling the vanishing 
of configurational entropy, and the presence of a thermo- 
dynamical Kautzmann transition. In Fig. 2, the solid and 
dashed lines represent the density evaluated analytically 
in the replica symmetric approximation, respectively in 
the liquid and solid (crystalline) phase, while the circles 
represent the density in the 1-step RSB ansatz. 

In the present paper, we investigate the dynamical be- 
havior of the model with z = 6, corresponding to the 
three dimensional case. We perform two types of Monte 
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FIG. 2; Density p as a function of f3p, calculated analytically 
(lines and circles) on the random graph with 2 = 6 (from 
Ref. [12]). Solid line represents the liquid phase, while the 
dashed line the crystalline phase. Circles: density in the 1- 
step RSB ansatz. Triangles: chemical potential as a function 
of the density, evaluated in the simulations performed in the 
present paper on a TV = 10^ random graph. 



Carlo dynamics on the system, canonical (fixed den- 
sity) and grand-canonical (variable density). The grand- 
canonical dynamics is given by the following algorithm. 

(1) Pick up a site at random. If the site is occupied by 
a particle, then pick up at random one of the z nearest 
positions (the other z — 1 positions inside the same site 
and the nearest position of the nearest site) and, if it is 
allowed by the hard core constraints, move the particle 
in the new position. 

(2) Pick up a site and a position at random. If the 
position is occupied by a particle, destroy the particle 
with probability exp{—p/T). If the position is empty, 
create a new particle if it is allowed by the hard core 
constraints. 

(3) Advance the time by 1/A'', where N is the number 

of sites. 

In the canonical dynamics step 2 is missing. Near the 
dynamical transition, where the probability to make a 
move is small, it is convenient to implement smart Monte 
Carlo algorithms like the "A/'-fold way" [14], that save 
computer time Iciaving the dynamics unchanged. 

We consider graphs composed by A'^ = 10^ sites (cells), 
each having z = 6 neighboring sites. First of all, we 
extract a regular random graph with fixed connectivity 
z. Starting from the empty graph, we perform grand- 
canonical dynamics with a slow cooling rate dp,/dt, un- 
til we reach the desired density. Therefore, we per- 
form canonical dynamics, skipping at least 10^ Monte 
Carlo steps for thermalization. All the quantities are 
averaged over 16 different realizations of the random 
graph. The relation between density and chemical po- 
tential can be evaluated in the canonical dynamics by 
Pp. — ln(Ap/ (Nf)), where Np is the number of particles, 
and (Nf) is the mean number of free positions, that is 
of positions where a particle might be created without 
violating geometric constraints. In Fig. 2, the correspon- 
dence between density and chemical potential in simula- 
tions is shown by the triangles. The values fall nicely on 
the curve corresponding to replica symmetry, indicating 
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that the system is in the equilibrium Uquid phase for all 
density considered. 




FIG. 3: a) Density-density relaxation function (j){t) as a func- 
tion of t/Ta, for densities (from left to right) p = 0.78, 0.80, 
0.81, 0.815, 0.817, 0.818, 0.819, 0.82. Inset: 4>{t) as a function 
of t. b) Fit of the function for the four densities nearest 
to the transition with the function defined in Eq. (7). 
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FIG. 5: a) Triangles: relaxation times Tq, defined as the time 

given by the fit e"*-*^^"' . Circles: time t* of the maximum 
of the dynamical susceptibility. Lines: fitting functions o|p — 
Pel"', with = 0.821 ± 0.001 and 7 = 2.75 ± 0.05. b) 
Maximum x(**) of the dynamical susceptibility as a function 
of the time t* . The straight line is a power law x(i*) ~ i*^ 
with exponent x = 0.36, equal within errors to I/7. 




FIG. 4: Dynamical susceptibility X'ii't) <^ function of time, 
for the same densities as Fig. 3a. Inset: The same plot in 
logarithmic scale, for the four densities nearest to the transi- 
tion. The straight lines are the power laws with exponents 
H — 0.29 and /x = 0.50, corresponding to exponents a and b 
of MOT. 

We then studied the dynamical functions. For each 
value of the density, we computed the quantity 



N 



P{z - P) 



(3) 



as a function of time, where n,i{t) = 0, 1 if site i is empty 
or occupied at time i, ai(t) = 0. . .5 is the position of 



the particle inside site i at time t, and p is the particle 
density. The density-density relaxation function of the 
system is defined as [(</>(t))], where (• • •) is the (equilib- 
rium) time average, while [• • • ] is the average over the 
graph configurations. The result is plotted in the In- 
set of Fig. 3a. We have fitted the long time tail of the 
functions with the Kohlrausch- Williams- Watts stretched 
exponential function /e~(*/'^") , and found that / and 
[3 are nearly constant for different densities, ranging re- 
spectively between 0.66 and 0.70, and between 0.74 and 
0.77. This shows that the time-temperature superposi- 
tion principle of MCT is well verified in the whole range 
of densities considered, as showed in the main frame of 
Fig. 3a. 

The dynamical susceptibility is defined as 



XA{t) = N[mf)-m)?]- 



(4) 



Fig. 4 shows the behavior of X4(*)- As expected, it has 

a maximum at a time i*, and both the time t* and the 
maximum Xi{i*) grow when approaching the transition. 
Note that the long time limit of Xi{i) is connected to 
a spin glass (static) susceptibility, lim(^ooX4(^) = XSG, 
where 



XSG 



(5) 
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and {n()niSa-„.oSai,a-) is the probability that there is a par- 
ticle in site and position (Tq = 0, and a particle in site 
i and position ct^ — a. As it is apparent from Fig. 4, the 
spin glass susceptibility remains small when approach- 
ing the transition, which means that static correlations 
between particles remain short ranged. 

The relaxation times Tc, extracted from the fit of 
{{(l){t))] for all the densities considered, and the time t* of 
the maximum of the dynamical susceptibility X4(^)i 
plotted in Fig. 5a respectively as triangles and circles. 
Both Ta and t* can be fitted with a power law, with the 
same (within errors) critical density p,. = 0.821 ± 0.001 
[15], and exponent 7 — 2.75 ± 0.05. Note however that 
the critical behavior is observed only very near to the 
transition, for \p — pc\ < 0.005, as anticipated in Ref. 
[11]. 

Prom the value of 7, we can extract the values of the 
exponents a and b that describe the approach and de- 
parture from the plateau and the exponent parameter A 
[7], 

^ 2a 2b' r(l-2o) r(l + 26) ^' 

The corresponding values arc a = 0.287 ± 0.004, b = 
0.50 ± 0.01, A = 0.785 ± 0.006. 

The MCT of the glass transition predicts that in the 
/3 regime the relaxation functions can be fitted with the 
function 

m = f + h^/\^\g{t/U), (7) 

where t„ = \a\~^ , a = p — pc, and ^(a;) is a specific func- 
tion depending on the exponent parameter A (for more 
details see Ref. [7]). We have fitted the relaxation func- 
tions with the function (7), for the four densities nearest 



to the transition, fixing a = p — pc with pc = 0.821, 
and extracting /, h and A from the fit (see Fig. 3b). 
The resulting parameters were / = 0.712 ± 0.003, h = 
0.95 ±0.15, A = 0.787 ±0.010. The parameter A is there- 
fore in agreement with the exponent 7 extracted from the 
power law divergence of the relaxation times Ta- 

We have then checked the predictions of MCT on the 
shape of X4(i)) and the divergence of Xii^*) ^ a function 
of t* [8, 9]. In the Inset of Fig. 4, we plot X4(t) on a loga- 
rithmic scale. As shown by the straight lines, during the 
early and late /3 regime (before it reaches the maximum), 
X4(i) can be well fitted by power laws t**, with p, = 0.29 
and p ~ 0.50, corresponding to the values of a and 6, as 
predicted by MCT. 

Finally, in Fig. 5b we show the dependence of the max- 
imum Xi{t*) as a function of the time t* . The values can 
be fitted with a power law Xi{i*) ^ i*^, with an expo- 
nent X = 0.36, compatible within errors with the value 
1/7, where 7 is the exponent of the divergence of times 
and t* . Also in this case, the critical behavior is observed 
only for \p - Pc\ < 0.005. 

In conclusion, the dynamical behaviour of the model 
on the random graph is consistent, near the transition, 
with the Mode Coupling Theory of the glass transition. 
In particular we have verified the quantitative prediction 
which relates the dynamical susceptibility with the decay 
of the density-density correlation function. Note that, 
due to the simplicity of the model, we have been able to 
consider a large system size. This is important because 
small systems suffer from large finite size effects near the 
transition, where the asymptotic behavior is expected. 
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